%% process_mesh
function mesh = process_mesh(pde, h, basis_type)
if h > (pde.domain(2)-pde.domain(1)) || h > (pde.domain(4)-pde.domain(3))
    error("h is too large.");
end
mesh.basis_type = basis_type;
mesh.Nx = (pde.domain(2)-pde.domain(1))/h;
mesh.Ny = (pde.domain(4)-pde.domain(3))/h;
[node, elem] = process_square(pde.domain, mesh.Nx, mesh.Ny);
[mesh.node, mesh.elem, mesh.edge, mesh.bde, mesh.bdn] = process_refine(node, elem, basis_type);
mesh.Np = size(node,1);
mesh.Ne = size(mesh.elem,1);
mesh.Ned = size(mesh.edge,1);
mesh.Nb = size(mesh.node,1);
mesh.Nlb = size(mesh.elem,2);
end

%% process_square
function [node, elem] = process_square(domain, Nx, Ny)
[X, Y] = meshgrid(linspace(domain(1), domain(2), Nx+1), linspace(domain(3), domain(4), Ny+1));
node = [X(:), Y(:)];
Np = size(node,1);
k = (1:Np-Ny-1)';
cut = (Ny+1)*(1:Nx);
k(cut) = [];
elem = [k, k+Ny+1, k+1;
    k+1, k+Ny+1, k+Ny+2];
elem = sortrows(elem);
end

%% process_refine
function [new_node, new_elem, edge, new_bde, new_bdn] = process_refine(node, elem, basis_type)
Np = size(node,1);
Ne = size(elem,1);
Ned = Np+Ne-1;
elem_ed = zeros(Ne,3);
edge = zeros(Ned,2);
edge_el = edge;
for e = 1:Ne
    n1 = elem(e, 1);
    n2 = elem(e, 2);
    n3 = elem(e, 3);
    tempe = [n1, n2;
        n2, n3;
        n3, n1];
    if e == 1
        edge(1:3,:) = tempe;
        count = 3;
        elem_ed(e,:) = [1, 2, 3];
        edge_el(1:3,1) = 1;
    else
        count1 = 0;
        for m = 1:3
            flag = 1;
            for i = 1:count
                if tempe(m,:) == [edge(i,2), edge(i,1)]
                    flag = 0;
                    elem_ed(e,m) = -i;
                    edge_el(i,2) = e;
                    break;
                end
            end
            if flag == 1
                count1 = count1 + 1;
                new_idx = count + count1;
                edge(new_idx,:) = tempe(m,:);
                elem_ed(e,m) = new_idx;
                edge_el(new_idx,1) = e;
            end
        end
        count = count + count1;
    end
end
bde = [edge_el(edge_el(:,2)==0,1), edge(edge_el(:,2)==0,:), find(edge_el(:,2) == 0)];
keyMap = containers.Map("KeyType", "double", "ValueType", "any");
for i = 1:size(bde,1)
    keyMap(bde(i,2)) = i;
end
new_bde = zeros(size(bde));
new_bde(1,:) = bde(1,:);
for i = 2:size(new_bde,1)
    idx = keyMap(new_bde(i-1,3));
    new_bde(i,:) = bde(idx,:);
end
bdn = unique(new_bde(:,2:3), "stable");
switch basis_type
    case "P1"
        new_node = node;
        new_elem = elem;
        new_bdn = bdn;
    case "P1b"
        center = zeros(Ne,2);
        for e = 1:Ne
            V = [node(elem(e,:),:)];
            center(e,:) = mean(V);
        end
        new_col = (Np+1:Np+Ne)';
        new_node = [node; center];
        new_elem = [elem, new_col];
        new_bdn = bdn;
    case "P2"
        midpoint = 0.5*[(node(edge(:,1),1) + node(edge(:,2),1)), (node(edge(:,1),2) + node(edge(:,2),2))];
        new_node = [node; midpoint];
        new_elem = [elem, abs(elem_ed)+Np];
        new_bdn = reshape([bdn'; (new_bde(:,4)+Np)'], [], 1);
    otherwise
        error("Invalid basis type.");
end
end